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SUMMARY 

Some  FORTRAN  programs  have  been  written  in  order  to  apply  the  Finite 
Element  Method  to  the  solution  for  low  Reynolds  number,  incompressible 
flows  around  a  Joukowski  aerofoil,  with  emphasis  on  the  generation  of  grids. 
These  programs  serve  as  evaluation  tools  and  as  a  first  step  in  a  planned 
longer-term  study  of  the  Finite  Element  Method  as  applied  to  fluid  flow 
problems. 
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1.  INTRODUCTION 

The  problems  of  external  flows  around  aircraft,  internal  flows 
inside  propulsion  units  and  also  related  problems  like  structural  designs 
all  require  more  and  more  accurate  solutions  while  incorporating  state  of 
the  art  features  involving  more  and  more  complex  geometries  and  loadings. 
This  is  steadily  becoming  beyond  the  reach  of  analytical  solution  methods 
and  thus  necessitates  the  use  of  new  computational  methods.  One  of  these 
is  the  Finite  Element  Method. 

The  Finite  Element  Method  was  initially  developed  and  used  by 
Zienkiewicz  'll]  for  elasticity  problems  and  is  at  the  height  of  its 
development  at  the  time  of  writing.  The  method  is  being  studied  theore¬ 
tically  as  well  as  being  applied  to  a  broader  and  broader  range  of  problems; 
its  applications  are  found  in  solid  mechanics,  fluid  mechanics,  electro¬ 
magnetics,  etc.. .A  reasonably  short  coverage  of  the  method  can  be  found  in 
the  book  by  Fenner  [3].  r 

The  FORTRAN  programs  given  at  the  end  of  this  Memo  have  been 
written  as  an  application  of  the  Finite  Element  Method  to  a  basic  fluid 
flow  problem,  namely  the  incompressible,  low  Reynolds  number  flaw  about  an 
aerofoil.  A  substantial  part  of  these  programs  deals  with  the  automatic 
generation  and  efficient  plotting  of  the  grids  surrounding  a  Joukowski 
aerofoil. 


The  work  reported  here  initiates  an  activity  which,  it  is  hoped, 
will  eventually  lead  to  a  capability  for  aerodynamic  estimation  beyond 
that  currently  possible  with  analytical  or  empirical  methods. 

2.  THEORETICAL  BACKGROUND 

The  Finite  Element  Method  is  a  method  applicable  to  linear 
problems,  initially  developed  for  elasticity  whereby  nodal  displacements 
in  a  continuum  are  solved  for  a  given  system  of  nodal  forces  provided 
these  nodal  forces  are  known  for  every  individual  mode  of  nodal 
displacement.  In  other  words  the  method  sets  up  a  large  system  of 
equations,  with  displacements  d^’s  as  unknowns,  forces  fj's  as  the  given 
constants  on  the  right  hand  side  of  the  equations,  and  the  coefficients 
a^j's  of  the  unknowns  are  written  down  from  the  knowledge  of  nodal  forces 
resulting  from  each  individual  nodal  displacement. 

Our  large  system  of  equations  thus  takes  the  form: 

*lldl  +  a12d2  + -  ”•  +  •ijAi  *  £1 

*21dl  +  *22d2  + -  *  •  •  +  *^**0  “  f 2 

(1) 


*nldl  +  *n2d2  +  "• 


"*  +  *1®^  *  £n 


2. 


where  a—  are  the  forces  generated  by  the  unit  displacement  dj  while  other 
displacements  are  kept  to  zero.  A  high  speed  computer  is  then  used  to 
solve  this  system  of  equations  when  the  f  j '  s  are  given.  The  above  system 
is  of  the  banded  type  due  to  the  process  of  forming  discrete  elements  to 
represent  the  continuum.  For  the  individual  displacement  of  each  node  the 
resulting  non-zero  nodal  forces  only  exist  at  its  ismediate  neighbour  nodes, 
each  of  which  must  share  at  least  one  comnon  element  with  the  node 
considered.  The  problem  of  boundary  conditions  where  displacements  instead 
of  forces  are  given  on  some  boundary  nodes  is  easily  dealt  with  by 
replacing  the  force  equations  for  these  nodes  with  the  simpler  displacement 
equations  which  have  the  known  values  of  displacements  on  the  right  hand 
sides  and  the  nodal  displacements  on  the  left  hand  sides.  A  special 
method  is  then  used  to  solve  this  banded  system:  considerable  saving  in 
computation  is  the  direct  result  of  the  purposely  designed  banded  character 
of  the  system. 

This  Finite  Element  Method  is  readily  adapted  to  the  problem  of 
Stokes  flow  where  the  above  displacements  are  replaced  by  fluid  velocities 
at  every  occurrence  (see  [3],  p.  16).  The  only  remaining  problem  is  the 
introduction  of  inertial  forces  as  a  kind  of  perturbation  to  the  viscous 
forces,  thus  retaining  the  linear  character  of  the  problem. 

The  central  problem  in  any  Finite  Element  Method  program  is  its 
efficiency  and  stability.  The  programs  given  at  the  end  of  this  Memo  are 
used  in  evaluating  these  two  aspects  for  a  certain  number  of  different 
schemes  available. 

A  substantial  part  of  the  programs  given  at  the  end  of  this  Memo 
deal  with  the  automatic  generation  of  a  grid  of  triangular  elements  around 
a  Joukowski  aerofoil.  The  grid  is  generated  by  wrapping  an  initial 
rectangular  mesh  consisting  of  triangular  elements  around  a  circular 
cylinder,  this  is  followed  by  a  conformal  transformation  to  modify  the 
cylinder  into  an  aerofoil  section.  Provisions  are  made  for  varying  the 
chord,  camber,  thickness  of  the  aerofoil  as  well  as  its  angle  of  attack. 

The  nodal  distribution  can  also  be  varied  to  put  more  nodes  in  the  boundary 
layer  than  in  the  far  field  uniform  flow. 

In  generating  the  mesh  the  numbering  of  the  nodes  affects  the 
bandwidth  of  the  system  of  equations  to  be  solved.  However  this  should 
not  alter  the  computation  speed  of  the  subroutines  for  solving  the 
equations.  The  grid  with  triangular  elements  can  surround  each  of  its 
nodes  by  the  least  number  of  six  adjacent  nodes,  hence  its  use  is  advanta¬ 
geous  in  reducing  the  bandwidth  of  the  system  of  equations. 

The  coordinate  tranformation  process  can  sometimss  invert  a 
finite  sized  element,  i.e.  change  its  oriented  boundary  into  the  opposite 
direction.  This  can  be  avoided  by  reducing  the  size  of  the  local  elements 
and  by  checking  the  sign  of  the  area  of  every  element  after  each  coordinate 
transformation . 
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DESCRIPTION  OF  THE  PROGRAMS 


Program  PEM  sets  up  the  grid  geometry  through  the  subroutines 
MESH3,  NDDFY3 ,  N0DFY4 ,  ELAREA  and  then  follows  the  standard  finite 
element  method  described  in  Section  2.  The  matrix  [a^j]  is  set  up  by 
the  subroutine  STIFF  with  the  boundary  conditions  catered  for  by  sub¬ 
routine  BCS.  The  system  of  equations  is  then  solved  by  subroutine  ELDON. 
The  solution  velocities  are  then  used  to  calculate  the  inertial  forces 
in  INERT  and  modify  the  nodal  values  of  forces  in  the  next  iteration  for 
a  more  accurate  solution  of  velocities.  Outputs  are  through  subroutine 
MSHOUT,  FEMDUT,  OUTPUT  and  SOLPLT. 

Some  features  of  these  programs  are: 

-  Since  the  program  FEM  originated  from  an  elasticity  problem 
it  is  necessary  to  set  the  Poisson  ratio  v  to  0.49  to  simulate  an 
incompressible  flow,  if  this  value  of  v  is  set  to  0.5  exactly  the  system 
has  a  number  of  infinite  coefficients  unless  the  problem  is  reformulated 
with  one  third  of  the  equations,  which  are  redundant,  removed.  Here  the 
approximate  value  v  -  0.49  is  used  to  avoid  the  complication  (see  [3], 
p.  154). 


-  Subroutines  STIFF  and  ELIKXN  are  for  fully  populated  matrices. 
They  are  used  here  only  for  quick  production  of  early  results  and  have 
been  replaced  in  subsequent  work  by  those  more  suitable  for  banded  matrices 
which  require  much  less  cocgwter  memory  and  time. 

-  A  seam  line  is  generated  at  the  trailing  line  of  the 
aerofoil  by  the  MESH,  M0DFY3,  M0DCT4  subroutines  to  put  a  number  of  nodes 
there. 


-  In  plotting  the  grid,  two  nodes  of  any  side  of  an  element 
are  joined  if  their  ordinal  number  increases  in  the  anticlockwise 
direction  around  the  element.  This  process  makes  the  plotting  computation 
linearly  proportional  to  the  nvmber  of  eleamnts  (hence  nodes)  and  avoids 
plotting  any  line  twice. 

-  Although  the  numbering  of  elements  and  nodes  does  not  effect 
the  speed  of  the  solution  routines  it  does  affect  the  speed  of  plotting 
the  grid. 

4.  RESULTS 

Figure  1  is  the  basic  rectangular  grid  consisting  of  triangular 
elements.  The  numbering  scheme  for  the  lTlelamants  and  105  nodes  is  self- 
evident.  This  whole  grid  is  then  wrapped  around  a  circle,  es  in  Figure 
2,  with  e  scaling  effect  to  put  more  nodes  nmer  to  the  inside.  The  grid 
than  undergoes  a  Joukowski  transformation  with  circulation  added  to 
became  the  grid  in  Figure  3.  The  Finite  Element  Method  ie  then  applied 
to  the  flow  field  using  this  grid.  The  resulting  velocity  field  for  a 
Reynolds  ntaber  of  1.2  (based  on  wing  chord)  is  plotted  in  Figure  4. 
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The  results  so  obtained  are  as  eapected  and  improvements  on 
the  method  are  under  study. 
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FI0.1  BASIC  RECTANGULAR  GRID  WITH  TRIANGULAR  ELEMENTS 


FI0.2  GRID  AROUND  A  CYLINDER 


FIG. 4  INCOMPRESSIBLE  FLOW  AT  REYNOLDS  NUMBER  1.2  AROUND  A  JCUKCWSKI  AEROFOIL 
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PROGRAM  FEN 


PROGRAM  FOR  FINITE  ELEMENT  ANALYSIS  OF  TUO- DIMENSIONAL  FLOW 
ABOUT  A  JOUKOUSKI  AEROFOIL,  USINE  CONSTANT  STRAIN  TR1AN6ULAR 
ELEMENTS. 

DIMENSION  XS0K300) 

REAL  NU 

COMHUN  /CAESH/  NEL,NNP,X<150),YU5v),AI<270),AJ<270>, 
AN(270>,BI(270),BJ<270) ,BN<270) ,AR£A<i70),NPI<270), 

NP J( 270 ) , NPK(  270 ) ,HOUT ,OLTY ( 270 ) 

/CMPAR/  NXPT,NYF  T  ,A,f> 

/CSTIF/  SSTIFH300, 300 ) ,nF i  i  <50 )  ,NPB2 l 50)  ,U(  1 50 ) ,0(150) , 
FX( 150) ,FT(150)  ,NBP1 ,NIP2 
/CMATL/  RO 

DATA  NROU,NCOL 7300,300/ 

L  GENERATE  MESH  DATA. 

TYPE  10 

10  FORMAT ('  INPUT  THE  NUMBERS  OF  POINTS  ALONG  THE  X  AND  Y  AXES/,/, 

RESPECTIVELY,  ALSO  THE  MESH  DATA  OUTPUT  CONTROL  PARAMETER  ,/, 

(0  IF  NO  MESH  DATA  TO  BE  PRINTED  OUT,  1  IF  ALL  MESH  DATA  TO  ',/, 
BE  PRINTED  0UT)s(3D  ) 

CALL  HESH3 
CALL  M0IFY3 
C  CALL  MODFYA 

C  COMPUTE  THE  ELEMENT  6E0METRIES. 

CALL  ELnREA 

C  OUTPUT  THE  MESH  BATA. 

CALL  MSHOUT 

C  PLOT  THE  MESH. 

CALL  nPL0T2 

C  SET  INITIAL  VALUES  01  STIFFNESSES,  EXTERNAL  FORCES,  AND 

C  VELOCITIES. 

BO  4  1*1 ,2*NNP 
BO  4  J*1 , 2*NNP 

4  SSTIFFTI , J)»0. 

BO  5  1*1, NNP 
FX< I )=0. 

FY( I )*0. 

U(I)*0. 

5  V(I)*0. 

C  ASSEMBLE  THE  6L0BAL  STIFFNESS  MATRIX. 

CALL  STIFF 

OPEN  (UNI  T  *?,*•!(.£* 'SS  TIFF.  TAP',  MODE* 'BINARY' ) 

UNITE  ( 9 )  Sb T IFF 
ILuitvONlT*? ) 

C  OUTPUT  THE  STIFFNESS  MATRIX. 

CALL  SHOUT 

C  INPUT  DENSITY,  »0. 

TYPE  46 

66  FORMAT < '  INPUT  THE  DENSITY,  RO:<F)> 

ACCEPT  47, kj 

67  FORMAT i F  i 


i 


OPEN (UN  IT *3, FILE* 'FEHOUT.DAT') 

UNITE ( 3,23) 

FORMAT ( '  SOLUTION  TO  FEM  PROGRAM' ) 

CL0SE(UNIT*3) 

SET  UP  ITERATIVE  SOLUTION  SO  LOOP. 

DO  71  L*1 ,3 

CALCULATE  INERTIAL  FORCE  PER  ELEMENT,  AND  DISTRIBUTE  AM0N6ST  THE 
ELEnttiTAL  NOSES. 

0P£N(UNIT*?,FILE*'SST1FF .TMP' .NODE*' BINARY  , 

REA b  (?)  SSTIFF 
CLOSE (UNIT*?) 

CALL  INERT 

APPLf  THE  BOUNDARY  CONDITIONS. 

CALL  ICS 

DO  »0  1*1 , NBP1 +NBP2 
IFd.LE.NBP1 )  N*NPBI  ( I ) 

IFd.bT.NBP1  )N=NPB2(  I-NBP1  ) 

FX(N)*U(N) 

FT(N)*V(N) 

DO  80  IC0L*1 ,2*NNP 
SSTIFF! 2*N-1 , ICOL ) *0 . 

SSTIFF(2*N,IC0L>*0. 

SSTIFF (2*N-1 ,2*N-1  )«1. 

SSTIFF ( 2*N, 2»N ) *1 . 

CONTINUE 
DO  21  1*1, NNP 

SSTIFF  <2*1-1 ,2*NNP+1 )*FX(I  > 

SSTIFF(2*I,2«NNP+1 )*FT(I) 

SOLVE  THE  SrSUH  OF  EQUATION*. 

MEQN*2*NNP 
TYPE  25, NEON 

FORMAT!  EN1ERIN6  ELlftlN* ,  NEON*' , 16) 

CALL  ELIMIN(SSTIFF,XS0L,HEQN,«R0U,NC0L,DET,AHEAN) 

OUTPUT  THE  SOLUTION. 

BO  36  1*1, NNP 
Ud)*XS0L(26l-1> 

V(I)»XS0L(2*I) 

CALL  FEMOUT 
CALL  OUTPUT 

CONTINUE 

OPEN! UNIT*?, FILE* 'SSTIFF.  THP') 

CLOSE (UNIT*?, DISPOSE»'DELETE'> 

PLOT  THE  VELOCITIES. 

CALL  SOLPLT 
END 

SUBROUTINE  MEsHj 

SUSPR06RAM  TO  READ  OR  6ENCRATE  A  MESH  OF  TRIANGULAR  FINITE 
ELEMENTS.  USED  IN  CONJUNCTION  UITN  H0Nfi3,  THIS  VERSION 
8ENERAIES  A  CIRCULAR  NESH  OF  MAINLY  ISOSCELES  TRIAN6ULAR 
ELEMENTS.  NODES  ARE  NUMBERED  IN  CONTINUOUS  RIN6S,  AND  ELEMENTS 
ARE  NUMBERED  IN  COMPLETE  ROUS,  I.E.  UPRIGHT  AMD  INVERTED 
TRIANGULAR  ELEMENTS  CONSIDERED  ALTERNATELY. 

COMMON  /CMESH/  NEL,NNP,X<150),Yd5*>,Al<270),AJ.270), 

AK(270) ,11(270) ,DJ(270) , DM (270), AREA! 270) ,NPI< 270, , 

NPJiI70) , NPK ( 270 ; , ROUT, SLTY( 270) 

/CMPAR/  NXPT.NYPT, A,l 


C  INPUT  THE  NUMBERS  OF  POINTS  AL0N6  THE  X  AND  V  AXES,  ALSO  THE 

C  HESH  DATA  OUTPUT  CONTROL  PARAMETER. 

ACCEPT  51,  NXPT ,NYPT ,HOUT 
31  FORHAT(3!> 

C  COHPUTc.  AND  TEST  THE  NUMBERS  OF  NODES  AND  ELEMENTS. 

HODNY*NOD(NYPT,2) 

IF(nODNT.EO.O)  NNP»NYPT*<2*NXPT-t 1/2 

If  i nODNY.EQ.1 )  NNP*<NYPT-1  )*<2*NXPT-1  1/2+NXPT-1 

NEL=<NYPT-1)*»2*NaPI-1) 

I F < NNP . LE . 1  SO . AND . NEL . LE . 270 )  60  TO  1 
TYPE  At,  NNP, NEL 

A1  FORMAT (30H0EXCESSIVE  SIZE  OF  MESH,  NNP  «,I3,8H,  NEL  », 15) 
STOP 

C  DEFINE  THE  NODAL  POINT  COORDINATES. 

1  1=0 

BO  3  IY*1 ,NYPT 
HODIY=HOD( IY,2) 

DO  2  IX*1 ,NXPT-1 
1*1  +  1 

X( I )*FL0AT(IX-1) /FLOAT (NXPT- 1) 

Yd)=FL0ATdY-1) /FLOAT (NYPT-1 ) 

2  IF(M0DIY.E0.0.AND.IX.6T.U  X(I)=X(I)-0.5/FLOAT(NXPT-1 ) 
IFfMODIY.Efl. 1  1  GO  TO  3 

1*1  +  1 

Y(I)*Y(I-1> 

X  < I ) - 1 . -0.5/FL0A1 (NXPT-1 1 

3  CONTINUE 

C  DEFINE  THE  NUMBERS  OF  THE  THREE  NODES  OF  EACH  ELEMENT. 

N=0 

NTEL«NYPT-1 
DO  4  IY*1,NYEL 
NXEL*N*PT-i 

IF(H0D(IY,2).EB.0>  NXEL*NXPT 

BO  4  IX*1 , NXEL 

N-M+1 

C  NNEU  REPLACES  H  AS  THE  HUMBER  OF  THE  ELEMENT 

HNEU«2*H 
NPI<HNEU)*M 
NPJ(HNEN)*NPI(MNEW>+t 
NPK ( MNEU)«NPI (NNEU 1+NXPT 
IF ( IX. Efi. NXEL ) HP J(MNEU) *NPJ< NNEU ) -NXEL 

4  IFdX.EO.NXPT  )NPK(HNEU)  *NPK(  NNEU) -NXPT  ♦  1 
M1*M 

DO  5  IT*1,NYEL 
NXEL-NXPT 

1F(H0B(!Y,2>.EQ.0)  NXEL*NXPT-1 
DO  5  IX*1 ,NXEL 

C  HNEU  REPLACES  M  hS  THE  NUMBER  OF  THE  ELEMENT 

H*M+1 

HHE«*(M'H1)*2-1 
NPI(MNEW>*(MNEU+1  )/2 
MP J ( NNEU ) *NP I (MNEU 1+NXPT 
NPR(HNEU)*NPI(MNEU)+NXPT-1 
IF(IX.E8.NXPT>MPI(MMEU)»NPUHHEU)-MXEL+1 
3  IF(IX.EO.NXEL )NPJ(MNEU/*NPJ.nNEV)-NXEL 

RE! URN 
END 


SUBROUTINE  M0DFT3 


SUBPROGRAM  TO  MODIFY  A  MESH  TO  SUIT  A  PARTICULAR  PROBLEM. 

THIS  VERSION  ADAPTS  A  SQUARE  MESH  10  A  CIRCULAR  hRC  (KING). 

COHHON  /CHESH/  NEL,NNP,X( 1 30)  ,Y( 150) , A I (  270 ) , AJ<270), 
hK<270),BI<2/(ij,BJ(270)iIM27G),AREA(270),NPI(270), 

NPJl  .’70)  ,NPK(2?0) ,  MOU  t , QLT  H  270 ) 

/CNPAR/  NXPT ,NYPT , h, ft 

INPUT  THE  SCALE  FACTOR  AND  THE  HESH  DIMENSIONS. 

TYPE  41 

FORMAT ('  INPUT  THE  SCALE  FACTOR  AND  THE  MESH  DIMENSIONS: !3F) ' ) 
ACCEPT  51,S,A,B 
FORMAT (JF ) 

TEST  FOR  ACCEPTABLE  BASIC  MESH. 

IF!M0D!NXPT,8).£Q. 1 .ANS.M0B(NYPT,2).£0. 1 )  60  TO  1 
TYPE  61 

FORMAT! '  HESH  UNSUITABLE  FOR  PRESENT  MODIFICATION') 

STOP 

PERFORM  FIRST  MODIFICATION  OF  i  COORDINATES. 

1*0 

HR*AL06(S)*(NYPT-1 ) 

DO  2  1*1 ,*NP 

C*Y  ( I ) «ALOG( 5) * ( NYPT- 1 ) 

Ifia.NE.1.0)C=lLXP(C)-1.0)/(EXP(HR)-1.0) 

IF ( S.EQ. 1 . 0 )C*f ( I ) 

Y!I)*C 

CONTINUE 

PERFORM  SECOND  MODIFICATION  TO  INTRUDUCE  CURVA1URE. 

PI*4.*ATAN< 1 . ) 

BO  3  1*1, MNP 
R*!A+(B-A)*Y! I )  )/2. 

PMI*X 1 1 )*2.0*PI 
X(I)*R»C0S(PHI)+0.5*B 
Y ( I )*-R*SlN!PHI )*0.3*B 


MODIFY  COORDINATES  OF  POINTS  NEXT  TO  THE  END  POINTS  OF  THE 
OUTERMOST  CIRCUMFERENTIAL  RON. 

I2*NNP-NXPT*3-!NXPT-U/4 

I1*NNP-NXPT*l 

BO  11  Ml, 4 

AR6-P1/4. -FLOAT <K)*PI/2. 

RM«aiRl(2.)*L0S(AR6> 

RX2*S0RT( 2. )*SIN( AR6 ) 

11*11 + (NXPT-1 ) /4 
I2*I2+!NXPT-1>/4 
I  TEMP* 1 1 
*1-12 
I2*ITEMP 

X( II )*B/2.»<1 .+RN1 ) 

T(I2)»B/2.*(1.+RK2) 

CONTINUE 

IF(NXPT.EQ.P)60T0  22 

II *NNF -NXPT+3- (NXPT-1 )/4 

I2*NNP-NXPT*1 

>0  22  X-1,4 

AR6*PI /4. -FLOAT! N)*PI/2. 

RXt*SSRT!2. >*C0S(ARB) 

RR2*$fiRT ( 2. ) *SIN i AA6 > 

1 i *  1 1 ♦ tNXPT-T )/4 
I2*I2*!NXPT-1)/4 


L  DEFINE  AnD  TEST  NEU  TOTAL  NUMBERS  OF  NODES  AND  ELEMENTS. 

I=NNF 

NNP»NNP* l NXPT-1 1/4+1-2 
M«NEL 

NEL«NEL+2*< (NXPT-1 >/4+1 )-4 
I F ( NNP .  LE .  1 50 .  AND .  NEL .  LE .  270  >  GO  TO  4 
TTPE  62,NNP,NEL 

62  FORMATS  EXCESSIVE  SIZE  OF  MESH,  NNP  NEL  *',I5> 

STOP 

C  IEFINE  THE  COORDINATES  OF  THE  ADDITIONAL  NODES. 

4  IXMAX*(NXPT-I T/4+1-3 
DO  A  IXa1 , IXNAX 
!■!♦! 

II»I1+IX 

IF(  IX.6T. ( (NXPT-1 )/4+1-3)/2i  60  TO  5 
X(I)*X(II> 

IF (MODTK, 2) . ED.  1 )  X(  I  )=B/2.»(»  .♦RM  > 

Y(I)*Y<II) 

IF (HCD(k,2) . E 0 . 0  >  Y(I>*B/2.*<1.+RK1) 

60  TO  6 

5  X(I)»X(II-1) 

IF ( MOOT  K , 2 ) . EQ . 0 )  XT  I )«D/2.*< 1+RK2) 

r<i)*=rdi-i  > 

IF (nuD(K,2) . EU . 1 )  Y ( I )-l/2.*( 1+RK2) 

6  CONTINUE 
Y(NNP>*B/2.*D/2.»RK2 
X(NNP)*B/2.+l/2.*RK1 

C  DEFINE  THE  NODES  OF  THE  ADDITIONAL  ELEMENTS. 

H1«H 

DO  7  IX* 1 , IXNAX 
N*N+1 

NPI (H )*I1 +H-M1 -1 
NPJ(H)*NPI(N)+1 

7  NPK(M)*NNP-( (NXPT-l )/4t|-2)+N-HT 

N2»N 

IXNAX*IXHAX-1 
DO  8  IX*1, IXNAX 
N«N+1 

NPI(N)»I1*N-N2 

NPJ(Ai*NNP-( (NXPT-1 1/4+1-2 JtH-H2+1 

8  NPK(M)*NPJ(H)-1 

NPI (NEL )*NNP-( (NXPT-1 )/4»/2 
NPJ(NEL )*NPI(NEL )*1 
NPK(Nc.k)*NNP 
22  CONTINUE 

C  CORRECT  COORDINATES  OF  NODES  ON  NESH  PERIMETER  FOR  SHALL 

C  DISCREPANCIES. 

DO  36  1*1, NNP 

DIFFX*X(I)-D 

DIFFY*Y(I)-I 

IF ( ADS < XT  I >  > . L T . 1 .E-4)XiI>*0. 

IF(ADS( Y( I ) ) ,LT. 1 ,E-4)Y(I )*0. 

IF(ADS(DIFFX) ,LT. I ,E-4)X(i)*D 
34  IF(ADSlDIFFY).LT.1.E-4)T(I>*D 


RETURN 

END 

SUBROUTINE  N0DFY4 


C 


SUBPROGRAM  TO  TKaNSFORM  A  CIRCULAR  ARC  TO  A  JOUKOVSKI  AEROFOIL 


COMMON  /CMESH/  MEL ,HNP,Xl 150) ,T 1 150) , A1 (270) , AJ<  270) , 

.  AK<2;0),&ii270),IJ<270),»K<270),AREA(270>,NPI<270>, 

.  NPJ<  270  ),NPX<  270 ) ,  MOUT , 0LTT<270 ) 

.  /CMPAR/  NXPT ,NYPT ,  A.I 

C  INPUT  THE  TUO  ECCENTRICITIES. 

1 TPE  41 

41  FORMAT ( '  INPUT  THE  IUO  ECCENTRITIESi <2 F)'> 

ACCEPT  51 ,X1 ,Y1 
51  FORMAT (2Fi 

X1*X1*A 
ri*ri*A 

C  INTROOUCE  TRANSLATION  OF  AXES. 

00  1  1*1 ,NNP 

C  ITPASS  TRANSFORMATION  OF  OUTER  PERIMETER  NODES. 

IF( (X( I) .EO.O. ) .OR. <X( I ) .EO.D.OR.  ( Y( I ) .E6.0. ) .OR.  < Y(I ) .EO.O) ) 
.  60  TO  1 

C  INCORPORATE  CIRCULaiION  CORRECTION  TO  SA1ISFY  KUTTA  TRAILIN6 

C  ED6E  CONDITION. 

XL*X ( 1 ) 

YL*T ( I ) 

PHI1*ASIN(Y1*2./A) 

R*SORT  l  ( XI - 0 . 5*B ) * ( XL-0. 5*1 ) ♦ 

.  <YL-0.5*D)*(YL-0.5*B>) 

PM1*AC0S( (XL-0.5*0)/R) 

IF  <  < YL-0.5*D) .LT .0. )PM1*-PH1 
PHI*FHI+PH1 1 *A/ (i»2. ) 

XL*R»CQS(PHI )+D/2. 

YL*R*SIN(PMI !+l/2. 

Af *SSRT( ( .5*A)**2-Y1*Y1 )-X1 

XI-(XL-.5*»-X1i/A1 

ETA»(YL-.5*i-Y1)/A1 

C  INTRODUCE  JOUKOMSKI  TRANSFORMATION. 

XISTA*XI*< 1 .*1 ./(XI4XI*ETA*ETA>) 

ETASTA-ETA*( 1 .-1 ./<XI*XI4ETA*ETA) ) 

C  REVERSE  INITIAL  TRANSLATION  OF  AXES. 

IF((X(I).NE.O.).ANS.(X(I).NE.») 

.  XU)*A1*XI$TA4.5*D+X1 

IF<(Y(I).NE.O.).ANS.(T(I).NE.O>> 

.  Y<I)»AI*ETASTA+.34l+Y1 

1  CONTINUE 

RETURN 

END 

SUDROUTINE  ELAREA 

i  SUDPROSRAM  TO  CALCULATE  THE  ELEMENT  AREAS  AND  DUALITIES. 

COMMON  /CME8N/  NEL,NNP,X(150),Y(150),AI<27«),AJ(270), 

.  AN  <  270 ) ,  D I  ( 270 ) ,  D J  <  270 ) . IN ( 27# ) , AREA ( 270 ) f  NPI ( 270 ) f 
.  NPJ<270) ,MPR<270) , MOOT ,OLTT 1270) 

/CMPAR/  MXPTfMYPI,A,D 

CALCULATE  THE  ELEMENT  AREAS  AMI  DUALITIES. 

DO  30  M*1 ,N£L 
I«MP1 (M) 

J»NPJ(M< 

KMPfcin) 

X1»X(J)-X<1) 

X2>X(K>-XU> 


C 


ruY»j(-r<i) 

Y2«Y(K)-Y<1) 

X3*X(K)-X<J> 

T3«V<K)-r ( J> 

$L2*S0RT(X2*X2+Y2*Y2) 

SL3*SBR1(X3*X3+V3*V3> 

AREA<H)«(X1*Y2-X2*Y1 )/2. 

BLT I' t  hs*AREA(M)/(SL1  +  SL2+SL3)**2.«12.*S0RT (3. ) 

IF (AREA(M) .6T.0. )  60  TO  30 
TYRE  32,11 

32  FORMAT ( '  ELEMENT  ',15,'  HA5  NE6AT1VE  AREA') 

STOP 

30  CONTINUE 

RETURN 

END 

SUBROUTINE  HSHOUT 

C  SUBPROGRAM  TO  URITE  OUT  THE  GEOMETRIC  DATA  FOR  THE  MESH. 

COMMON  /CHESH/  NEC ,NNF,X< 1 50) , u t 50) , AI (270) , AJ(270> , 

.  AK<270),BI<270),BJ(270),BK(27G>,AR£A(270),NPli^/u), 

.  NFJ(270i ,NPK(270) , ROUT, BLTY( 270) 

.  /CMPAR /  N*PT,NfPT,A,B 

IF(NOUT.EO.O)  RETURN 

C  OUTPUT  THE  NUMBER  OF  ELEMENTS,  NODAL  POINTS  AND  COORDINATES. 

TYPE  M,  M£L,NMP,(I,X<1) ,Y<1) ,I*T  ,MHP> 

01  FORMAT ( '  6E0HETRIC  BATA  FOR  THE  MESH',//, 

.10X, 'NUMBER  OF  ELEMENTS  *',I4,//, 

.10X,  NUMBER  OF  NODAL  POINTS  «',I4,//, 

.'  NODAL  POINT  COORDINATES',//, 

I',4X,'X',iX,'Y',/ 

. (IX,I5,2F$.4>) 

C  OUTPUT  THE  ELEMENT  NODE  NUMBERS  AND  AREAS. 

TYPE  62,  (M,NP1 (n> ,ftPJ<H) ,MFKiH) ,AREA(H),ILTY(H) ,H«1  ,NEL) 

62  FORMAT (//,'  ELEMENT  NODE  NUMBERS, AREAS  AMI  OVALITIES  ,//, 

.  M',4X,'1',4X,'J',4X,'I('  ,6X,'AREA',7X, 'IUALITY' ,/, 

.<1X,4I5,1X,Et2.4,4X,F4.3>> 

RETURN 
END 

SUBROUTINE  HPL0T2 

C  SUBPR06RAH  TO  PLOT  THE  NESH. 

COMMON  /CHESH/  NEL,MNP,Xn50),Y(150),AI<27B),AJ(270), 

.  AK(270),BI(270),BJ<270) , BK ( 270 ), AREA ( 270 ),NP1( 270) , 

.  NPJ<270) ,HPK( 270) , ROUT, SLTY< 270) 

/CNPAR/  NXPT,NYPT,A,B 

•1NENSIDN  NFN7(4),XN(4),im><< 

0FEN(UNIT*1 ,FILE«'HSHPLT.0UT' ) 

C  kEPOSIUOM  PLn. 

CALL  PLOT ( 1 ,0. ,5. ,2) 

SCALE* 10. 


BO  61  A«I,NEL 
NPNT< 1 )*NPI <N) 


MPNT<2>«NPJ<N) 

NPNT(31*NPK<M) 

NPNT(4>«NPI<N) 

DO  71  K*l ,3 
I1*NPN1<K) 

12=HP»(T  <K*1  > 

XN(1>-X(I1> 

XN(2)«X(I2) 

XM(3>«0. 

XM<4)*1 ./SCALE 

YN<1)*YUU 

YN<2)*Y<I2) 

YR<3>-0. 

TN(4)*1 ./SCALE 

CHECK  UMETHER  PAIR  OF  POINTS  ARE  TO  IE  JOINED. 

IFLA6-0 

IF  < 1 1 .IT. 12)  IFLhu-1 

CHECK  iiriETHER  PAIR  OF  POINTS  LIES  ON  NtSH  OUTER  PtUnElEK. 
IF(((XN(l).E0.0. ) .AND. (aN(2).EQ.0.j).QR. 

<  (XN<  t ) .EO.D ) .AND. (XN(2).EO.D)).OR. 

<  <  YN(  1 ) .EQ.O. ) .AND. < YN<2> .CB.O. ) 1 .OR. 
((YN(1).EB.B).ANw.(YN(2).EU.D)))  IF  LAS* 1 

CHECK  UHETHER  PAIR  OF  POINTS  LIES  OR  HESH  INNER  PERIMETER. 
IF < (II .LE.NXPT-1 ) . AND. < I2.LE.NXPT-1 ) )  IFLA6*1 

JOIN  APPROPRIATE  POINTS. 

1FUFLA6.E0.1)  CALL  LINEd ,XN, TN,2, 1 ,0,46) 

CONTINUE 

1ETERH1NE  ELEMENT  CENTROIL  COORDINATES. 

I*NPI<N> 

J»MPJ<NJ 

K*NPK(M) 

AfclC^t  lX>I)*X(JMX(K)  )/3.  )*SChLL 
YELC-i »i<I)*Y(J)*Y<K) )/3. ) ‘SCALE 

DETERMINE  NOHDER  OF  DI6ITS  IN  M. 

NDN*INT< AL0010(FL0AT(M) I >+1 

ADJUST  ELEMENT  NUiilCR  COORDINATES. 

AXElC'XELC- . 5*FL0AT(NBM+1 )«.06 
AY£LC*YELC-.5*.07 

MURDER  ELEMENT. 

CALL  MURDER ( 1 ,AXELC,AYELC, .07  ,FLOAT<M> ,0. ,-1 ) 

XST AR*AXELC+FLOAT (MBH+1 )* .06- .03 

XST AR«XELC+FLOAT ( RDM ) * . 04 - . 03 

CALL  SYMBOL < 1 ,XSTAR,YELC,. 07,42,0. ,-i) 

CONTINUE 

MURDER  NODES. 

DO  SI  1*1  ,NNP 
Ta«X\I>*SCALE 
TY*T( I )«SCALE 

CALL  hUMDER« l.TX.TI,. 07, FLOAT* I ),0.,-1> 

CONTINUE 

CLOSE ( UNIT* 1 , FILE* 'ASHPIT .  OUT' ) 

RETURN 

END 


•UDR0UT1NE  STIFF 

SUDPROSRAM  TO  FORM  INDIVIDUAL  ELEMENT  STIFFNESS  MATRICES,  AND 


c 


ASSEMBLE  THE  OVERALL  STRUCTURE  STIFFNESS  MATRIX 


REAL  NU 

IIMENSION  IJK(3>,IMAT(3,A),0<3,3),ESTIFF<M> 

COMMON  /CMESH7  MEL ,NNP,X<  ISO) , T<  1 30) , AIi'270)  ,AJ(270> , 

.  AK (270), BI(270),IJ(270), IK (270), AREA (270), NF 1(270), 

.  NF  ji 270) ,NPK( 270) ,H0UT,0LTY (270) 

.  /CMPAR/  NXPT , NYPT , A ,  I 

/ CST I F /  SSTIFF( 300,300) ,MPI1 (50) ,NPI2(50 ) ,U( ISO) , V( ISO ) , 

.  FX(I50),FT( t50),MIP1 ,NkP2 

C  INPUT  THE  HAIEKIAL  PROPERTIES  OF  THE  ELEMENTS. 

TYPE  27 

27  FORMAT)'  INPUT  THE  VALUE  OF  E  AND  NU:(2F>') 

ACCEPT  37,E,Nu 
s7  FORMAT <2F) 

DO  SO  I*1,2*NNF 
DO  SO  J*t,2»NNP 
50  SSTIFF(1,J)«0. 

C  SET  UP  OVERALL  ASSEHILY  LOOP. 

DO  20  n-l ,NEL 

C  COMPUTE  THE  ELEMENT  6E0HETRIES. 

l-NPI(H) 

J»NPJ(M> 

K-NPK(H) 

AI(M)**X(J)+X(K) 

AJlM)«-X(K)+l(I) 

AK<M)«-X<I)*X<J) 

II<M)«Y< J)-Y(K) 

IJ<H)»Y<K)-Y(I) 

IK<H)"Y<I)-Y<J> 

C  STORE  THL  ELEMENT  NODE  NUMBERS  IN  ORDER  IN  ARRAY  I  JR. 

IJK(D-NPHM) 

IJK(2)«NPJ(M) 

IJK<3)«NPK(N)  ♦ 

C  FORM  THE  ELEMENT  DIMENSION  MATRIX,  IHAT . 

•0  7  IRE* 1 ,2 
BO  7  1CE-1.A 

7  DMAT(IRE,ICE)*0. 

IMAT(1,t )aBI(M> 

BMAT ( t ,3)*IJ(H) 

IMATO  ,5)*DK(h) 

IMAT(2,2)*AI(M) 

»MAT(2,4)«AJ(Mi 
DNA)(2,*)«AK(N) 

•0  8  ICE* l ,  A 

If (M0BUCE,2).E0.0)  8HAT<3,ICE>B»MAT<1,ICE-n 

8  IF (HOD (ICE,2).EQ.1>  DMAT(3,ICE)-»MAT(2,ICE+1 ) 

C  FORM  THE  ELASTIC  PROPERTY  MATRIX,  D. 

10  9  IRE*1 ,3 
•0  f  ICE*1 ,3 

9  t( IRE, ICE>>0. 

FACT»E/(( 1 .♦NU)*<1 .*2.*NU) ) 

B<I,T)«FACT»<I.-NU) 

»(1,2)»FACT«MU 
B(2, t )>)(1 ,2) 

D(2,2)*B(),T > 
l(3,3)«FACT*.5*(1.-2.*HU) 


C  FORM  THE  ELEMENT  STIFFNESS  MATRIX,  ESTIFF. 

S3  10  1=1,6 
10  10  J* 1 » 6 
ESTIFFI I, Ji*0. 

BO  10  1*1,3 
10  10  1<*1 ,3 

10  ESTIFF <  I,  J)«E$T IFF (l,J)«IMT(L,ll*Ba,K>*BHAT<K,J>*.2S/AREA<H> 

C  CONSTRUCT  OVERAU  STRUCTURE  STIFFNESS  MATRIX,  SSTIFF. 

10  2V  IR£=1 ,3 
BO  20  IC£*1 ,3 
IROU*I JK( IRE) 

ICOL*I JK( ICE) 

SSTIFF (2*IR0U'1 ,2*IC0L-1 )*SSTIFF(2»IR0U-1 ,2*IC0L- 1 ) 

.  ♦£STIFF(2*IRE-1,2*ICE-n 

SSTIFF(2*IR0o-1 ,2«JC0L )*SSTIFF ( 2*IR0U- 1 ,2* I  COL) 

.  ♦ESTIFF<2*IRE-1 ,2*ICE) 

SSTIFF  (2«IR£m*,2»1C0L-1  )  =  SSTIFF  t2*IRQU,2*IC0L-1 ) 

.  ♦ESTIFF i 2*1  RE, 2* ICE- 1 ) 

SSTIFF (2»IR0U,2«IC0U*SSTIFr  (2*IROU,2*ICOU 
.  ♦ESTIFF(2*IRE,2»ICE ) 

20  CONTINUE 

TYPE  5000 

5000  FORMAT (  STIFF  EXECUTES') 

RETURN 

END 

SUBROUTINE  STF0U1 

C  SUBPROGRAM  TO  OUTPUT  STRUCTURE  STIFFNESS  MATRIX. 

COMMON  /CHESH/  NEL ,NNP,X< 150) , Y < 1 50 ) , AX (270 > , AJ»270 ) , 

.  AK 1270) ,11 ( 270),  BJ(  270),  BK(  270) ,  AREA  1 270),  NH  (270/ , 

.  NPJ(270),NPK(270),M0UT,QLTY<2?0) 

.  /CUPAR/  ti*FT,NTPT,A,B 

.  ,CSUF/  SSTIFF (300, 304)  ,NPBt  (SO) ,NPB2<50) ,U( 150) ,V(  150), 

.  FX( 1 50 ) , FT ( 150) ,NBP1 ,NBP2 

TYPE  SI 

61  FORMAT ('  INPUT  THE  6L0BAL  STIFFNESS  NATRIX  OUTPUT  CONTROL',/, 

.  '  PARAMETER  (I  IF  THE  EHTIRE  STIFFHESS  MATRIX  10  BE  TYPED  ',/, 

.  '  OUT,  0  IF  HOME  OF  THE  STIFFHESS  MATRIX  TO  BE  TYPEB  OUT)  XI)  > 

ACCEPT  71, ROUT 
71  FORMAT ( I) 

IFOCOUT.EO.O)  RETURN 

TYPE  50 

50  FORHAT »/  /  , '  GLOBAL  STIFFNESS  MATRIX') 

BO  t  1*1 , 2*NMP 
TYPE  15,1 

15  FORMAT (/,  kuU  ,i*i 

IiFE  1 w, t SSTIFF <1,J),J*1 ,2*NHP) 

10  FORMAT  IGF 12.7; 

i  CONI i HUE 

RETURN 

end 

subroutine  inert 

c  subprogram  to  calculate  the  inertial  force  per  element,  and 

C  DISTRIBUTE  IT  AH0N6ST  THE  ELEMENTAL  NODES. 

CONNOR  /CHESH/  NEL,NNP,Xn30),Y(lS0),AI(270),AJ(270), 

.  AR(270),B1<27G>,BJ(270),IR<270>,AREA(270),NPI(276>, 

.  NFwi//0),NPR(270),NOUT,IUT(270) 

.  /CSTIf /  SSTIFf v 300,300) ,NPI1 (SO) ,NP|2(50> ,U( ISO) ,V( 150) , 


.  FX( 150),FY( 150),NBP1 ,NBP2 
.  /'Ctifc TL/  RO 

DIMENSION  NPNH4) 

ll U  50  1*1  ,HHp 

F*a>*o. 

so  rr<V‘0. 

DO  41  h*i ,NEL 
NPNTU>*NPl(H> 

NPNT(2)»NPJ<N> 

NPNT (3>*NPK(H) 

NPNT(4)*NPI(N) 

FEX*0. 

FEY*0. 

Du  81  K«t ,3 
I1«NPNT(K) 

I2*NPNT(K+1 ) 

D2«U( I2)-U(  II ) 

D3*V(11) 

D4*V(I2)-V(I1 ) 

RL1*X( 12)-X(11 > 

RL2*Y(I2)-YU1> 

Rl«S»ta  v RL1*Rl1»RL2*Kl2> 

RNX*RL2/RL 

RNY—RL1/RL 

FEX*FEX*RO*RL*(fii*(  (D1+D2)*RNX-*»ii*.5*D4)*RNi  >»D2*ii2*R((X/3. 

.  ♦( .5*D3«BA/3, )*RNY 1 1 

Ft  i-tLf*kO*liL»<B3*(  (B1  +  .5*D2)*RNX+ (D3+D4)*RNY )*B*» 1 1 . 5*D1+D2/3. ) 
.  *RNX+B4«RNY/3. ) ) 

81  CONTINUE 

DO  31  K*1,3 
I1«NPNT(K) 

FX(I1)*FX(I1)-FEX/3. 

31  FVill )*FY< 11 ) -FEY/3. 

41  CONTINUE 

RETURN 

END 

SUBROUTINE  DCS 

C  SUBPR06RAM  TO  APPLY  THE  BOUNDARY  CONDITIONS. 

CONNON  YCNCSK/  NEl,NNP,XUSOMt150>,AIt2?0>,AJl270,,. 

.  AX<270),II<270),IJ<270)  ,BK( 270), AREA) 270), NPK 270), 

.  NPJ<270 I  ,Nf'k  \  »  ,MflUT ,BLT i 

.  /CNPAR/  NXPT,NYFI ,A,t 

.  /CSUf/  SST IFF (300, 300)  ,MPB1  <5<it  ,NPB2<30i  ,ov IS® ) ,V( 150) , 

.  F  X ( 1 50 ) , F  Y ( 1 30 ) , NBP 1 , NBP2 

C  STORE  NUMBERS  OF  NODES  ON  INNER  NESH  BOUNDARY  IN  ARRAY  NP»1,  AND 

C  »ET  CORRESPONDING  VELOCITIES  TO  ZERO. 

1*0 

BO  10  N*1,NXPT-1 
1*141 
NPBUD-N 
ttt»)*0. 

10  V(N)*0. 

NBP1-I 

C  STORE  NUMBERS  OF  NODES  ON  OUTER  NESH  DbUNBAK)  IN  ARRAY  NPD2,  AND 

C  SET  CORRESPOND INS  VELOCITIES  TO  UNITY. 

1*0 


SO  20  N=1 ,NNf 

IF<<X<N).N£.0.).AND.(X<N).Nt.&).AND.rTiN/.NE.O.>.f*ND.d<N>.NE.B) 

.  )  bit  10  20 

1*1+1 
MPB2<  i )  =** 

U(N)=1 . 

V<N)*0. 

20  CONTINUE 
NBP2*I 

RETURN 

END 

SUIROUTINE  ELIMIM1 A, X,MEON,NROU,NCOL,DET  ,AHEAN) 

C  SUBPROGRAM  FOR  SOLVING  SIMULTANEOUS  LINEAR  EQUATIONS  BY  GAUSS  I AM 

C  ELIMINATION  U1TH  PARTIAL  PIV0TIN6. 

DIMENSION  A(NROU,NCOL ) ,X(NROU) 

DOUBLE  PRECISION  SUM 

NEQN*HEQN 

IF(NE8N.LE.NR0U.AND.NESN.LE.NL0L'  i )  bit  TO  t 
WRITE (6, 61 ) 

61  FORMAT (33H0ST0P  -  DIMENSION  ERROR  IN  EL  In I N ) 

STOP 

C  FIND  MEAN  COEFFICIENT  HA6N1TU0E. 

1  A«EAN*0. 

DO  2  1*1 , NEON 
DO  2  J*1 , NEON 

2  AMEAN*AME AN+ADS ( A ( I , J ) ) 

AHEAN*ANEAN/FLOAT (NEQN*NEOM) 

C  COMMENCE  ELIMINATION  PROCESS. 

JMAX«NE0N*1 
NE0NH1-NE0N-1 
DO  6  IEON*l , NEBNN1 

C  SEARCH  LEADING  COLUMN  OF  THE  COEFFICIENT  MATRIX  FROM  THE 

C  DIAGONAL  DONNVARDS  FOR  THE  LARGEST  ELEMENT  AND  RAKE  THIS  THE 

C  PIVOTAL  ELEMENT. 

IHIN*IE6N* i 

!MAX*!£8N 

BO  J  1*1MIN,NEUN 

S  IF  kAiai A( 1 , IE8N) ) .GT. ABS( A( IMAX, IEON) ) )  IMAX-I 

IF < IMAX.Efl. IEON)  60  TO  S 
DO  4  J*IEQNf JHAX 
AA*A< ICON, J) 

AdEON,  J)*AdNAX,J) 

4  AdMAX,  J)*AA 

C  ELIMINATE  *<IE0N>  FROM  EOUATIONS  dEONO)  TO  NEON,  FIRST  TESTIN6 

C  FOR  NONZERO  PIVOTAL  ELEMENT. 

5  IF (ABS( AdEON, IEGN)/AHEAN) .LT .  1  .E-S)  60  TO  10 

DO  6  2*IMIN,NE0N 

F ACT -At  1 , 1 £0M  > / A ( IEOH, IEQN) 

DO  6  J*IMIN, JHAX 

A  Ad,J>«Ad,J)  -FAC  T»A  dEON,  J ) 

C  SOLVE  THE  UPPER* TOIANOUL At  SET  OF  EOUATIONS  DT  BACK 

C  SUBSTITUTION. 

IF(ADS(A(NEON,NEON»/AMEAN).Li. i.E-01  GO  TO  10 
X ( NEON  >  * A  <  NEON , JMAX ) / A ( NEON , NEON  > 

DO  6  L*2,NE0N 
I*NEONd-t 


SUH-A ( I , JMAX  > 

DO  7  J=IP1 ,NEQM 
SUM*SUH-A(I,J)*X(J> 

X<  I ) =SUM/A( I ,  I ) 

EVALUATE  DETERMINANT  OF  COEFFICIENT  MATRIX  AND  COMPARE  UITH 
ORIGINAL  COEFFICIENTS. 

BETA*1 . 

DO  9  1*1 ,NEQN 
DETA*D£TA*A( I , I ) 

IET-DETA 

TTPE  72 

FORMAT ( '  ELIMIN  EXECUTED  ) 

RETURN 

DET*0. 

TYPE  66 

FORhAU'  ELIMIN  NOT  EXECUTED  ') 

REIOhft 

END 

SUDROU1 1NE  FEMOUT 

SUBPROGRAM  TO  STORE  THE  SOLUTION  BATA  IN  OUTPUT  DATA  FILE. 

COMMON  /CMESH/  «EL,NNP,X( 1 50 ) ,Y< 1 50) ,A1 ( 270) , AJ(270> . 

.  AK  <  270 ) , BI < 270 ) ,BJ(270),BK< 270) , AREA (270) ,NPI (270) , 

.  NPJ(270) ,NPK(270) ,HQUT,0LTY(270) 

.  /CMPAR/  NXPT,NYPT 

/CSIIF/  SST IFF (300,300) ,NPB1 (50) ,NPB2 (50) ,U( 150 ) ,V( 150 ) , 
.  FX( 1 50 ) ,FY< I 50/ ,NBP1 ,NBP2 


OPEN ( UNI T*3, FILE* 'FEMOUT. BAT', ACCESS* 'APPEND' ) 

WRITE (3, 11),(I,FX(I),FY(I),U(I),V(I>,I*I ,NNP) 

FORMAT (//,'  NODAL  FORCES  AND  VELOCITIES',//, 

.  '  N',6X,'FX',10X,'FY',11X,'U',11X,'V',/, 

.  <  X , IS , 4E 12.4) ) 

CL0SE(UNII*3, FILE* 'FEMOUT. DAT') 

RETURN 

ENB 

SUBROUTINE  OUTPUT 

SUBPROGRAM  TO  OUTPUT  THE  SOLUTION  DATA. 

COHnON  /CMESH/  MEl,NMP,X(130),Y(150),AI(270),AJ(270), 

.  AM270),BI(270),BJ(270),IK(270),AREA(270),NPI(270), 

.  NP J( 270 ), NPK (270), MOUT, OLT i(270 I 
.  /CMPAR/  NXPT,MYPT 

.  /CSTIF/  3»TIFF(300,300) ,NPI1 (50),NPB2(30),U( 150) ,V( 150) , 

.  FX( 1  SO ) ,F  Y ( 1 50 ) ,NDP1 ,NBP2 


TYPE  It , (I,FX(I),FY(I),U(I),V<I),I*1,NNP) 
FORMAT*//,'  NODAL  FORCES  AND  VELOCITIES',//, 
.  '  N',6X,'FX',10X,'FY',11X,'U',t!X,'V',/, 

.  (X,I5,4E12. A) ) 

RETURN 

ENB 


SUBROUTINE  SOLPLT 


SUB*R06RAM  TO  PLOT  THE  NODAL  VELOCITIES. 

CGHftoN  /CNESH/  NEL ,NNP, X( 1 50) , Y( ISO) ,  AI (270 ) ,AJ< 270 ) , 

.  AK(27G),SI(270),BJ(270)  ,BKi270>,AREAi27«/,NP 1(270 ) , 

.  NPJ(270),NPK(270) , MQUT ,QLT I i2/w) 

.  /CMPAR/  NXPT,NYPT,A,B 

.  /CST IF /  SST IFF (300,100) ,NPBl (50) ,NPB2( 50) , U ( 1 50 ) , V( 1 30) , 

.  FX<15Q),FY<150),NBP1,NDP2 

DIMENSION  XVECT (4),YV£CT(4),XN(4),TN(4) 

OPEN (UNIT-2, FILE«'SOLPLT. OUT') 

REPOSITION  PEN. 

CALL  PLOT (2,0. ,5. ,2) 

SCALE-10. 

00  f  1=1, NNP 

COMPUTE  VtLOLi f Y  VECTOR  COORDINATES. 

XvtCl ( I  )»X( I ) 

XVECT<2)=X(I)+U<I)/20. 

XVECT  <3)=0. 

XVECT  <  4  J  =  1  ,/SLAlE 
YVtLT ( 1 )*Y( I ) 

YVECT(2) *Y ( I )+V( I ) /20. 

YVECTv3)*0. 

YVECT (4)*1 ./SCALE 

PLOT  VECTORS. 

CALL  LINE12, XVECT, YVECT, 2, 1,0,46) 

XNODE=XVECT( 1)*SCALE 
YNODE=YVECT ( 1 )*SCALE 

CALL  SYMBOL (2,XN0SE, YNODt , .07 ,4,0 .  ,-1 ) 

CONTINUE 

PLOT  AEROFOIL. 

DO  51  1*1 ,NXPT-2 
XN( t )*X( I ) 

XN<2)»X< 1*1 ) 

XN(3)*0. 

XN(A)«1. /SCALE 
YN( 1  )»Y( I ) 

TN(2)=Y(It1 > 

YN(3)*0. 

YN<4)*1. /SCALE 

CALL  LINE(2,XN,YM,i, i,0,4o) 

CONTINUE 

XN( I )=X(NXPT-1 > 

XN(2)*X< 1 ) 

YN<1)«Y<NXFT-i> 

r«(2)=Y<1) 

CALL  LINE(2,XN,YN,2, f ,0,46) 

CLOSE < UN  I T«2,FILE*'$0LPLT. OUT') 


RETURN 

END 
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